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ABSTRACT 


~P  This  report  deals  with  the  problem  of  pitch  and  heave  damping 
of  a catamaran  with  small-waterplane-area-twin-hull  (SWATH)  con- 
figuration. A computer  program  has  been  developed  to  compute  the 
pitch  and  heave  damping  coefficients  of  SWATH  including  forward 
speed  effects  based  on  the  thin  ship  theory  developed  by  Newman 
(1959).  Calculations  of  the  damping  coefficients  for  several 
Froude  numbers  are  compared  with  Lee's  (1974)  experiments.  The 
results  show  that  damping  is  greatly  affected  by  the  Brard  para- 
meter, wV/g  when  the  forward  speed  effects  are  considered.  The 
effect  of  hull  distance  variation  on  damping  is  also  considered. 
The  results  show  that  damping  is  an  oscillatory  function  of  hull 
distance.  The  oscillatory  phenomenon  dies  out  as  the  hull  dis- 
tance increases,  and  also  as  the  Brard  parameter  increases.  . 
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I. 


INTRODUCTION 


Due  to  the  advantages  of  large  deck  area  and  stable  motion  charac- 
teristics in  heavy  seas,  small-waterplane-area-twin-hull  (SWATH)  type 
of  catamaran  has  become  an  interesting  subject  recently. 

In  recent  work  on  the  prediction  of  catamaran  motions  in  waves, 

Lee  (Pien  and  Lee  [11] 1 ) extended  the  strip-theory  to  catamaran  configu- 
rations. A fundamental  assumption  in  the  strip-theory  is  that  the  hydro- 
dynamic  characteristics  of  a ship  can  be  inferred  from  the  two-dimensional 
stripwise  characteristics  of  each  section  along  the  length.  The  effects 
of  forward  speed  are  included  in  Lee's  work  by  the  strip-theory  synthesis 
rather  than  by  a rigorous  introduction  of  forward  speed  effects  into  the 
hydrodynamic  boundary  conditions. 

The  results  based  on  the  strip-theory  are  generally  in  good  agreement 
with  the  experiments,  except  for  a pronounced  resonance  effect  in  the 
theory  at  a critical  frequency.  At  the  corresponding  wavelength  the 
motions  are  substantially  overpredicted  compared  with  experiments.  This 
resonance  seems  to  be  a consequence  of  the  presence  of  near-zero  damping 
at  zero  forward  speed  and  the  use  of  the  zero-speed  hydrodynamic  coef- 
ficients in  a strip-theory  manner.  This  problem  occurs  only  for  hull 
forms  where  bulbous  cylindrical  sections,  having  a small  waterplane 
area  and  a large  submerged  volume,  are  dominant. 

To  be  more  specific,  let  us  first  consider  a two-dimensional  thin-body 
section  with  a shape  y “ +h(z)  . Then  it  is  known  that  the  damping 

'Numbers  in  brackets  designate  Reference  at  end  of  paper. 
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coefficient  is  proportional  to  the  square  of  the  integral 


Kz  dh 


where  T is  the  draft  and  K is  the  wavenumber.  (See  Section  2.1).  It 
seems  apparent  that  for  a bulbous  form  where  dh/dz  changes  sign,  the 
above  integral  will  vanish  for  a suitable  combination  of  the  wavenumber 
and  hull  shape.  This  phenomenon  has  been  investigated  earlier  by  Motora 
and  Koyama  [7].  They  did  experimental  work  on  two-dimensional  wave- 
excitationless  forms,  which  have  similar  forms  to  the  SWATH  demi-hull. 
Frank  [4]  did  some  computations  on  heave  damping  for  various  bulbous 
cylinders.  In  both  studies,  it  was  shown  that  there  is  a critical 
frequency  at  which  the  damping  coefficients  vanish. 

In  three  dimensions  with  zero  forward  speed,  the  situation  is  changed 
because  the  damping  coefficient  depends  on  an  integral,  with  respect  to 
the  waveangle  0 , of  the  square  of  the  surface  integral 

/ / eKz  + iKx  CO80  3htx»z>  dz  dx 

J J H 7 


where  L is  half  length  of  the  body.  Since  the  damping  integral  is 
positive-def inite  for  all  0 , zero-damping  seems  impossible.  However, 
it  can  be  anticipated  that  there  might  occur  near-zero  damping  for  long 
cylindrical  vessels  having  essentially  the  appropriate  two-dimensional 
form  because  there  will  be  a dominant  waveangle  ^ as  suggested  by  strip 
theory  or  the  stationary  phase  approximation  of  the  three-dimensional 
damping  integral  for  KL  >>  1. 


Finally,  if  forward  speed  effects  are  included,  the  possibility  of 


zero  damping  is  even  less.  In  this  case,  the  three-dimensional  damping 
coefficient  is  again  proportional  to  the  square  of  a surface  integral 
similar  to  that  shown  above,  but  now  the  wavenumber  K is  no  longer  a 
single  constant,  but,  depending  on  the  value  of  the  Brard  number  a)V/g  , 
takes  on  either  two  or  four  discrete  values,  each  of  which  depends  on 
9 (See  Newman  [8]).  Thus,  since  the  square  of  the  surface  integral  is 
Integrated  over  a continuous  spectrum  of  K,  the  probability  of  zero 
damping  is  greatly  reduced. 

Calculation  of  the  three-dimensional  damping  coefficients,  including 
forward  speed  effects  without  the  assumption  of  strip  theory  is  so  compli- 
cated that  some  alternative  simplifying  assumptions  are  required.  Newman 
[8]  presented  such  a theory  for  "thin"  mono-hulls,  including  illustrative 
calculations  for  a mathematical  hull  form  for  which  experimental  data 
had  been  obtained  by  Golovato.  Subsequently,  Gerrltsma,  Kerwln  and  Newman 
[5]  presented  comparisons  of  the  same  theory  with  experiments  for  a Series 
60  hull  form.  In  both  cases  the  comparison  was  qualitatively  useful. 

In  this  report, Newman's  thin-ship  theory  is  applied  to  the  catamaran 
hulls,  especially  with  SWATH  configuration,  in  an  attempt  to  avoid  the 
deficiency  of  the  strip  theory  at  nonzero  forward  speed.  The  fundamental 
assumption  of  Newman's  thin-ship  theory  may  be  valid  for  these  hull  forms, 
since  they  are  thin  in  the  Important  region  near  the  free  surface.  Moreover, 
to  extend  the  practical  validity  of  the  present  results,  the  submerged 
cylindrical  hull  portion  will  be  treated  by  a modified  thin-ship  approach 
as  described  in  Section  2.3.  Based  on  Newman's  work  [8],  pitch  and  heave 
damping  coefficients  are  calculated  for  the  NSRDC  Model  MODCAT.  The  results 
are  compared  with  Lee's  [6]  experiments  for  various  Froude  numbers. 


IL.  DAMPING  OF  TWO-DLMF.NSIONAL  BULBOUS  CYLINDERS 


2 . 1 Damping  of  a Two-Dimensional  Thin  Body 

Consider  a thin  vertical  body  which  is  in  an  oscillatory  heave 
motion  on  the  free  surface  with  velocity  £ = V cos  Lit  . The  fluid  is 
assumed  ro  be  inviscid,  irrotational  and  of  infinite  depth.  Assume  the 
body  is  symmetrical  about  the  z-axis,  z is  positive  upwards  and  its 
hull  function  is  given  by  y = + h(z)  for  -T  < z < 0 , where  T is 

the  draft  of  the  body.  Then  the  velocity  potential  which  satisfies  the 

the  linearized  free  surface  condition,  in  the  region  of  positive  y , is 
known  to  be  [10]: 

0 

<p  = 2eKz  V sin(Ky-wt)  / — ^ eK^  d £ 

-T 

GO  0 

2 ,,  _ r r dh  -ky 

- — v cos  ut  / / -fir-  e 

0 -T 

(k  cos  kz  + K sin  kz)  (k  cos  k£  + K sin  kQ  ^ 

k(k2  + K2) 

(2. 

where  the  first  term  on  the  right  hand  side  is  related  to  the  outgoing 
waves  and  the  second  term,  to  the  local  disturbance. 

From  Bernoulli's  equation,  the  total  hydrodynamic  force  acting  on 
the  body  in  z-direction  is  obtained  as  follows: 

0 

Z - 2 / p cos(n,z)  ds 
-T 


10  - 


= 4oop  V cos  wt 


-T 


dh  Kz  , 

— : e dz 

dz 


4 t 

- cup  V sin  (at  / 

v 1 


1 2 


/ — (k  cos  kz  + K sin  kz)dz 


k^  + k^) 


dk 


= b33  V cos  wt  - a33  Vwsin  tat 

(2.2) 

where  b33  is  the  damping  coefficient  and  a33  is  the  added  mass  coef- 
ficient. From  the  relation  (2.2)  we  iiave  a formula  for  the  damping  coef- 
ficient : 


'33 


4<ap 


r _Kz  dh 

I 6 IT  dz 


(2.3) 


2 . 2 Extension  of  Thin  Ship  Result  to  Bulbous  Cylinder 

As  a means  to  check  the  validity  of  thin  ship  theory  to  find  the 
damping  coefficients  of  SWATH  configurations,  a two-dimensional  model 
(Fig.l)  is  considered.  In  this  section,  damping  of  a heaving  circular 
cylinder  with  a thin  vertical  strut  Is  examined  by  using  (2.3).  Define 
a polar  coordinate  system  (r,0)  with  the  origin  at  the  center  of  the 
cylinder,  0 being  positive  clockwise  starting  from  z-axis.  Let  h(z) 
be  the  hull  function,  rQ  the  radius  of  the  cylinder,  and  0O  the  angle 
between  z-axis  and  the  bottom  of  the  strut.  If  dh/dz  = 0 for  z > -T 
+ r0(l  + cos80)  , the  integral  in  (2.3)  becomes: 
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where  I 


r Kz  U n 
I e --- — dz  = 
' dz 


t Kz  dh 
/ e 


-T 

0. 


dz 


dz 


= / 


Kr.  cos  0 - K(T-rfl ) 


0 r0  cos0  d0 


TT 


-KT0  r Kr0  cos0  Q 

- r0  e 0 J e u cos0  d0 

0, 


0C 


r0  e KT°  {irl^Krp)  - j eKr°  cos®  cos©  d 0 } 


j (z)  is  the  modified  Bessel  function  defined  by 


1 j z COS0 
0 


1,(2)  - -±-  / e“  cosO  d0  . 


i 


Ceometry  o*  Bulbous  Cylinder 


"ISP 


Assuming  Kr0  sin0o  <<  1 , the  integral  in  (2.4)  becomes: 
0o 

/ eKr°COS0  cos0  d0  - 60  e*^0  + O(0o2)  . 

o 

Hence,  the  integral  in  (2.3)  is  reduced  to  be: 


0 

/ 

-T 


Kz  dh 
dz 


dz 


% 


TTr 


e"KTo  I j (Kr ) + e e-K(To-ro  ^ 


(2.5) 


In  a sense,  the  first  term  on  the  right  hand  side  of  (2.5)  represents 
the  effect  of  the  circular  cylinder  and  the  second  term,  the  effect  of 
the  thin  strut.  As  an  extreme  case  where  e is  negligibly  small  and 
Kr0  -*•  0 , then  Ij(Kr0)  ^ y Kr„  and  (2.5)  is  approximated  by: 
o 


Kz  dh 


dz 


dz 


<v  - 


TTr 


K e'KTo 


(2.6) 


-T 


which  represents  only  the  effect  of  the  submerged  cylinder.  Then  the  far 
field  behavior  of  the  velocity  potential  from  (2.1),  using  (2.6),  becomes: 


4>  ^ 2V  eKz  sin  (Ky  - cut)  (-  ^ TTr02  K e"KT°)  (2.7) 

Based  on  (2.5),  the  damping  coefficients  are  calculated  for  several 
cylindrical  forms.  The  results  are  compared  with  Frank's  [4]  in  Fig.  2 
and  with  Lee's  [6]  in  Fig.  3.  Agreements  are  good  for  only  low  frequency 
range.  Poor  agreements  in  the  frequency  range  of  practical  importance 
indicate  that  we  need  to  modify  the  thin  ship  results  for  the  circular 
cylindrical  part  which  is  actually  not  thin  and  for  which  the  thin  ship 
assumption  may  break  down. 
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2 . 3 Correction  for  the  Submerged  Cylindrical  Part 

In  order  to  account  for  the  discrepancies  in  damping  coefficients 
resulting  from  the  thin  ship  approximation,  Eqn.  (2.5)  will  be  modified 
for  the  effect  of  the  submerged  cylindrical  part  by  comparing  the  far 
field  behavior  of  the  velocity  potential  (2.7)  with  that  of  the  known 
solution  of  the  circular  cylinder.  Consider  a circular  cylinder  of  radius 
rQ  located  at  z = -T„  which  is  oscillating  with  velocity  C = Vcosoot  . 
Assuming  T0  » r0  and  using  the  same  coordinate  system  as  in  Fig.  1, 
the  velocity  potential  is  known  to  be: 

<p  - -V  rQ2  coswt  C°^  ^ (2.8) 

It  is  known  (e.g.  (13.31)  in  Wehausen  & Laitone  [13])  that  the  velocity 
potential  of  a two-dimensional  source  log  r located  at  (0,-0 

behaves  at  infinity  like  QeK^z+^+iy^  where  Q is  the  source  strength 
and  K is  the  wave  number.  The  corresponding  behavior  of  the  potential 

of  a vertical  dipole  — ^ will  be  QKeK^Z+^+iy^  since  the  velocity 

potential  of  a vertical  dipole  can  be  obtained  by  differentiating  the  source 
potential  with  respect  to  C • Then  it  can  be  readily  shown  from  (2.8)  that: 


Comparing  (2.9)  with  the  limiting  case  of  the  thin  ship  result  (2.7),  they 
differ  by  a factor  of  2 . Thus  it  is  reasonable  to  modify  the  thin  ship 
results  (2.5)  by  doubling  the  effect  of  the  circular  cylindrical  part: 


f Kz  dh  . nt  „ 

I e — - — dz  ■ — 2tt 

• n •* 


Note  that  the  laat  term  in  unmodified  reflecting  the  importance  of  the 
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waterplane  area  in  the  low  frequency  limit.  The  same  arguments  will  also 
apply  to  the  three-dimensional  case  when  we  perform  the  integration  of 
the  hull  function  which  has  the  similar  form  to  that  in  (2.5).  Damping 
coefficients  based  on  the  modified  thin  ship  result  (2.10)  are  plotted 
in  Fig.  2 and  Fig.  3. 


III.  EXTENSION  OF  THE  THIN-SHIP  THEORY  TO  CATAMARAN 


3.1  The  Interaction  Effects  Between  the  Two  Hulls 


Consider  now  three-dimensional  catamaran  hulls  with  cross  sections 
having  essentially  the  same  as  two-dimensional  bulbous  cylinders 
separated  by  a distance  2b  described  in  Fig. 4.  Newman  [8]  developed  a 
theory  of  the  damping  of  a thin  ship  by  the  analysis  of  energy  radiation 
in  surface  waves.  After  an  asymptotic  expansion  of  the  Green's  function, 
pitch  and  heave  damping  coefficients  of  a single  hull,  by  separation  of 
the  energy  components,  were  found  to  be: 


3C0VB2  />-  jPi2  + (tK-1) 4sgn(TK-l) 

17  _oo  1 P 2 2 + Q2  [ (tK-1)  4 - K2]^ 


(3.1) 


where  is  the  pitch  damping  coefficients  and  Zw  , the  heave  damping 

coefficients  in  the  linearized  equations  of  motion.  Based  on  (3.1)  the 
interaction  effects  between  the  two  hulls  of  the  catamaran  will  be  included. 


In  order  to  simplify  the  problem,  only  a first  approximation  to  the  hull 
interaction  effects  in  considered;  i.e.  the  source  distributions  of  each 
separate  hull  are  linearly  superposed  as  in  the  wave-resistance  theory 
for  catamarans  [3].  Then  the  expansion  of  Green's  function,  equation  (66) 
in  [8],  is  slightly  changed  by  the  effects  of  the  twin  hulls 
in  the  form  of  exp  (-iA^  b sin  u)  + exp  (iA^  b sin  u) } . Then  the  inte- 
gration of  the  hull  function  (P.j.Q.j,)  » equation  (69)  in  [8],  is  multi- 
plied by  the  same  factor  for  catamaran  hulls: 

(P,,Q.)  = (P.,Q.)  4 i u n * 2 cos (A . b sin  u)  (3.2) 

i xi  catamaran  i,xi  single  hull  i 
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J . 2 Jamp lng  .>t  SWATH  Catamaran 

Us  1 rig  the  same  notation  as  [8],  the  damping  coefficients  of  pitch 
and  heave  in  the  linearized  equations  of  motion  for  catamaran  hull 
incorporated  with  (3.2)  will  be: 


. « . 2P^J;  ['  . 

fh’+Q.1] 

r*  2 , /n  2 

• (^--.l.)‘,^n(jK-lj_  |4cos2[vb/(TK_1)4  _k2] 

TT  _oo 

(p2  + Q2  J 

[ (xK-1) 4 - K2]2  ( 

where  M * pitching  moment 
Z * heave  force 
q =*  pitch  angular  velocity 
W * heave  velocity 
p = fluid  density 

to  = circular  frequency  of  oscillations 
T = odV/g 

V = w2/g 

8 = beam-length  ratio 

b = half  hull  separation  distance 


‘l-//  ( 

* s i 


3h  r 5h 

C 35  * 3C 


Sitl(vK£)  ev^('rK_1) 

cos 


dK  d c 


dK 
(3.3) 


(3.4) 


0. 


^ = I!  -It-  !*"(vK£)  eVi;(TK~1)-  d£  dc 


cos 


(3.5) 


The  prime  in  equation  (3.3)  denotes  that  only  the  intervals  of 


(iK-l)”*  - KJ 


K,  < K < K 
1 - - 2 


are  to  be  included  in  the  integration.  The  intervals 


and  K3  < K < K4  are  omitted,  where 


20 


K 


l 


K 


2 


> 


[ (2t-1)  + 


(1 


— - [ (2t+1)  + (1+4t)^] 


(3.6) 


(3.7) 


Thus  the  integral  in  equation  (3.3)  can  be  decomposed  in  the  following 
forms  according  to  the  values  of  T : 

C 

Kj 

/ F(K)  dK  - I if  l“0 

K2  (3.8) 


/'  F (K)  dK  = J 

-OO 


K3  » 

/ + / +/ 

— ( oo  kf  1/ 

2 ; 


F (K)  dK 


I2  + !3  + \ 


if  0<x<y 

(3.9) 


F(K)  dK  = I5  + 16  if  T>\ 

(3.10) 

where  F(K)  denotes  the  integrand  of  (3.3).  Numerical  procedure  to 
evaluate  the  integrals  in  equations  (3.8),  (3.9)  and  (3.10)  are  described 
in  the  subsequent  chapter. 


/ +/ 


** 


3 . 3 The  Effects  of  Hull  Distance  on  Damping 

It  can  be  easily  anticipated  that  the  damping  of  a catamaran  is 
affected  by  the  hull  separation  distance  as  well  as  the  forward  speed  due 
to  the  interactions  between  the  generated  waves  and  hulls.  The  equation 
(3.3)  shows  that  damping  is  an  oscillatory  function  of  the  hull  distance. 
For  a simple  case  of  heave  damping  with  zero  speed,  the  asymptotic 
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behavior  of  the  damping  when  b is  large  can  It  idily  derived. 


From  (3.3),  the  heave  damping  integral  .hen  T=0 


Z = - 
v 


— f (p.2  + Q,  2)  2[cos(2vb\i'i -K1  I ->  1]  dK 

1 " 2 {l-K 

(3.11) 


In  order  to  examine  the  asymptotic  behavic:  ) for  large  b , 

consider  only  the  oscillatory  integral  involving  the  cosine  term: 

I = ±_2a -l&H  cosCZvb^2)  dK  (3.12) 

1 1 fl-K2 

From  (3.5),  (P22  - Q 2)  can  be  shown  to  be  an  even  function  of  K 

Then  (3.12)  becomes: 


- 1 f p ^ ■+•  0 ^ ^ f — TT 

1=2'  . Jl  JL-  co8(2.Vb\l-Kz)  dK. 

1 \jrF-P 


(3.13) 


3y  a change  of  the  variable,  x “Vl-K2  : 

lt  = 2 , jL  ^ ^2  F.  cos(2vbx)  dx 


..2r  \ (p22  + q22) 
c — 


f l+x 


cos  (2vbx) 


dx 


fi- 


cos(2vbx) 

2 : j.  x . 


dx. 


yl-x 


(3.14) 


where  f(x)  is  a regular  function  between  (0,1)  By  a successive 
integration  by  parts  (see  Copson  [2]),  it  can  be  shown  that  the 
leading  order  asymptotic  behavior  as  b -*■  +°°  will  be: 


li  % cos(2vb  " 


(3.15) 


IL 
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IV.  NUMERICAL  PROCEDURE 


4 . 1 Zero  Forward  Speed  Case 

When  T=0  , K2  and  K3  become  -1  and  +1  respectively.  Thus 
equation  (3.3)  is  reduced  to  the  following  simple  form: 


2pqjvg^ 


/ dK 


V + Q:2 


V + Q2‘ 


+ Q 2 


■(  1_-  {4cos2  [vbVl-KZ  ]} 
1 -R2  L ■> 


(4.1) 


The  integral  in  (4.1)  is  easily  evaluated  using  the  Gauss-Chebyshev 


quadrature  formula: 


/ -^=-  dK  = l w±  F(K±)  + En 
-l  V 1-K2  i=i 


(4.2) 


where  Kj  are  the  roots  of  the  mth-degree  Chebyshev  polynomial,  so  that 


K_^  = cos 


! wi  = F/m  ; and  En  is  an  error  term. 


Then  (4.2)  is  simplified  to: 


-!  Vl-K2 


dK  = — — y F 

TT1  - “ 


2(1-1)* 


(4.3) 


Damping  coefficients  are  non-dimensionalized  by  the  quantity  pV\/g/L  for 
heave  and  by  pVLV gL  for  pitch  where  V is  a displaced  volume.  Based  on 
(4.1)  and  (4.3),  pitch  and  heave  damping  coefficients  for  a single  and 
twin  hulls  are  calculated  for  the  NSRDC  model  MODCAT . Results  are  plotted 
in  Fig. 6 and  Fig. 7. 
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4 . 2 Nonzero  Forward  Speed  Case 

If  the  forward  speed  effects  are  included,  not  only  the  finite 
integral  but  also  the  semi-infinite  integrals  should  be  evaluated.  The 
semi-infinite  integrals  are  quite  involved  due  to  the  highly  oscillatory 
cosine  term.  Consider  the  two  cases  separately  according  to  the  values 
of  T . 

(a)  0 < T < 1/4 

In  this  case,  the  damping  integral  is  decomposed  in  three  different 
ranges  as  in  (3.9).  The  finite  integral  I3  can  be  treated  in  the  same 
manner  as  the  zero  speed  case  by  an  appropriate  change  of  the  variable  of 
the  integration.  Rewriting  the  finite  integral: 


= / (pi2  + Qi2) 


(tK-1)  sgn(lK-l) 
V(TK-1)“  - K2 


4cos2  [vbV(xK-l)1*-  K3]  dK 


/ j(Pt2  + Qt2) 


t2V(k4-k)(k-k1) 


: [vb^F-l)4-  K2]}~ — -- — — - 

J V(K-Ky)(K3-K) 


Denoting  the  expression  in  the  braces  above  by  F (K)  , 


k3  F,(K) 

I3  - / dK 

V(K-K  )(K3-K) 

By  a linear  change  of  the  variable  of  integration, 


(4.4) 


K3-K2 


(K-K.)  + 1 


(4.5)  reduces  to: 


*3  " / 


k3-k2 


k3+k2\ 


\IT~^ 


(4  . b) 


L 
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where  we  can  use  Gauss-Chebyshev  quadrature  formula. 


The  semi-infinite  integrals  in  (3.9)  and  (3.10)  can  be  treated 
essentially  in  the  same  way.  In  order  to  facilitate  the  integrals,  we 
manipulate  in  the  following  way: 


I = / (P.2  + Q.2)  • 4cos2 [vbV (tK-1) **-  K2] 

2 L,  * 1 fCrK-l)“  - K2 


dK 


= / (P  2 + Q,  2)  CLL-lj.NgnCtK-l^  .2[cos(2vb\fcnr-T ) + l]-^_ 

T2V(K4-K)(K3-K)(K2-K)  \^-k 


K,  r 

2 / J(P  2 + Q 2)  (JK-!)  Bgn_([K-l)  _ , cos[2\W(TK-l)4-  K2] 


dK 


l/(K4-K)(K3-K)(K2-K) 


/KrK 


+ 2 / (P  2 + Q 2)  (T-K~1)-  

^ _ 2 ll/f,  ir\  /,r  „\ 


dK 


T 2 \[0VK)  (K3-K)  (k2-k)  V Kj  -k 


- 2(1  + I.)  (4.7) 

a b 

where  I denotes  the  first  integral  with  an  oscillatory  integrand  in  (4.7) 
a 

and  , the  second  one.  Concerning  only  the  first  integral, 


-/ 


(Pj2  + ql2) 


(TK-l)‘,sgn(TK-l) 


T2i/(K„-K)(K3-K)(K2-K) 


— cos 


2vb ^(iK-1) K2] 


dK 


\f^ 

(4.8) 


Dominant  contributions  to  the  integral  come  from  the  vicinity  of  Kj  since 
the  oscillations  get  faster  as  j K ) increases,  thus  cancelling  out  effectively. 
Equation  (4.8)  can  be  written  in  the  simple  form: 


F2(K) 

\T VK 


dK 


(4.9) 
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where  F2(K)  denotes  the  expression  in  the  braces  in  (4.8).  Then  by  a 
change  of  the  variable,  x = ^TCj-K  : 


1=2  / F, (K.  - x2)  dx 


(4.10) 


The  semi-infinite  integral  is  subdivided  into  an  infinite  number  of  finite 


ones  of  an  interval  A , so  that: 


oo  oo  A(n+i ) 

1=2/  F2(Kj-x2)  dx  = 2 l / F2(K1-x2)  dx 


(4.11) 


n=o  An 


However,  due  to  the  oscillatory  nature  of  the  integrand  in  (4.8),  we  should 
be  very  careful  in  choosing  the  size  of  the  interval  A . In  order  to 
account  for  the  change  in  period  of  an  oscillation,  the  interval  A is 
chosen  as  the  period  of  cos  (2vb\](xK-l)1*-  K^)  . Then  the  sub-integrals 
in  (4.11)  are  performed  over  these  periods.  As  n increases,  the  contri- 
bution of  the  sub-integral  gets  smaller.  Thus  the  integration  can  be 
performed  to  a desired  accuracy  by  controlling  the  upper  limit  of  n . 

After  the  value  of  A is  determined,  the  sub-integral  1^  in  (4.11)  for 
a given  value  of  n , becomes: 


\ = / F.O^-x2) 


(4.12) 


where  a and  b are  the  lower  and  upper  limits  of  the  sub-interval. 
Then  by  a change  of  the  variable  of  integration  to: 


2x  - (a+b) 
b - a 


(4.13) 


Equation  (4.12)  reduces  to: 


b - a 

2 


/ F2i  Ki  " 


(4.14) 
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We  now  have  an  appropriate  integral  form  for  which  we  can  use  the  Gauss- 


Legend  re  quadrature  formula: 

1 n 

/ F(z)  dz  = l w.FCz^ 
’ i=o 


-1 


(A. 15) 


where  w^  are  weight  factors  given  by: 


l n 

= / IT 

-1  j = 0 

j^i 


zi  - zj 


dz 


(A. 16) 


and  z^  are  the  roots  of  the  Legendre  Polynomial  Pn+1  (z)  . The  roots 
z i and  the  weight  factors  w-j  for  several  values  of  n are  listed  in 
Table  1. 

For  the  integral  I4  , the  same  procedure  can  be  used  as  in  I2 
except  for  a different  change  of  the  variable  of  the  integration.  From 
(3.9), 


h = / MK> 


dK 


(A. 17) 


K-K„ 


where 


F4(K)  = (P±2  + Q12) i 


- { 4cos  2 rvb  ^(-tK:d-_  kt1 

\j( K-KJ  (K-K2)  (K-Kj)  1 1 L 


By  changing  of  the  variable,  x = \J  K - K : 


I4  = 2 / F4 (K4+x2)  dx 


(A. 18) 


Analogous  to  (A. 11),  we  subdivide  the  semi-infinite  integral  into  finite  ones 


00  A (n+i ) 

I4  = 2 l / F4(K4-Hc2)  dx 


(A. 19) 


n=o 


An 


By  the  same  change  of  the  variable  (A. 13),  each  sub-integral  in  (A. 19)  re- 
duces to: 
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. ! 


h - Fo K- + 


z(b-a)  + (a+b) 
2 


dz  (4.20) 


where  the  Gauss-Legendre  quadrature  formula  (4.15)  and  (4.16)  may  be  used, 
(b)  T > 1/4 

When  T > 1/4  , we  have  two  integrals  Is  and  I6  similar  to  12 
and  I4  . briefly  repeating  the  same  procedure  as  in  I2  and  I4  : 

K, 


and 


where 


h ! F5<K) 


h = / 

k4 


dK 


\[k3-  k 


dK 


fK- 


(4.21) 


(4.22) 


(TK-l)^gn(TK-l)  1 f. 


F 5 (K)  - (P.2  + Q.2)  -I  l 4cos2(vb^^y^ 

1 1 t\|[(tK-1)2+  K](K4-K)  [ [ L 


(4.23) 


and 


Fh (K)  = (Pj2  + Qi  ) 


(tK-1)  ‘*sgn(TK-l) 
iffOlK-l)2!-  K]  (K-K  ) 


4cos 2 j Vb  ^(tK-1) ** - K: 


(4.24) 

By  changing  the  variables,  x = yK3  - iT  for  (4.21)  and  x » ^TR  - K4  tor 
(4.22),  (4.21 ) and  (4.22)  reduce  to: 


and 


I 5 = 2 / F5(K3-x2)  dx 
0 


h - 2 / F6(K4-x2)  dx 
o 


(4.25) 


(4.26 


where  we  can  use  the  Gauss-Legendre  quadrature  formula  for  each  sub-integral 
after  dividing  the  semi-infinite  integrals  into  finite  ones  as  done  for 


I2  and  I4  . 
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F{z)dz  £ n-,f(s,) 


Roots  (z,) 


Weight  Factors  («■,) 


Two-Point  Formula 

n = 1 


±0.577.35  02691 

89626 

Three-Point  For  inula 
n = 2 

1.00000 

ooooo 

0 (XXX) 

0.00000  00000 

ooooo 

0.88X88 

88X88 

88889 

±0.77459  66692 

41483 

Four-Point  / in  inula 
n =3 

0.55555 

55555 

55556 

±0.33998  10435 

84856 

0.65214 

5)548 

62546 

±0.86113  63115 

94053 

Fire-point  Formula 
n 4 

0.34785 

48451 

37454 

0.00000  00000 

ooooo 

0 56X88 

88X88 

888X9 

±0.53846  93101 

056X3 

0.47862 

86704 

99.366 

± 0.90617  98459 

3X664 

Six- Point  1 in  mula 

..  C 

0.23692 

68850 

56189 

oc 

cc 

ri 

o 

-H 

83197 

n 7 

0.46791 

39.145 

72691 

± 0.66120  93X64 

66265 

0.36076 

15730 

48139 

±0.93246  95142 

03152 

Ten-Point  Finmiilu 
n 9 

0.17132 

44923 

79170 

±0.14887  43389 

X 1 6.3 1 

0.29552 

42247 

1475.3 

±0.43339  53941 

29247 

0.26926 

67 1 93 

09996 

±0.67940  95682 

99024 

0.21908 

63625 

15982 

±0.86506  33666 

88985 

0.14945 

13491 

50581 

±0.97390  65285 

17172 

Fifteen-Point  Formula 
n — 14 

0.06667 

13443 

08688 

0.00000  00000 

OOOOO 

0.20257 

82419 

25561 

±0.20119  40939 

974  35 

0.19843 

1485.3 

2711 1 

: 0.39415  13470 

7756.3 

0.18616 

KKKII 

15562 

± 0.57097  21726 

08539 

0. 1 667.6 

92058 

1 6994 

± 0.72441  77.313 

60170 

0 1 3957 

06779 

26154 

± 0.84820  65834 

10427 

0.10715 

92204 

67 1 72 

± 0.93727  3.392.4 

00706 

0.07036 

60474 

88108 

1 0.98799  25180 

204X5 

0.03075 

32419 

96117 

. Roots  of 

the  Legendre  Polynomials 

Pn+!(Z) 

and 

the 

Factors  for 

the  Gauss-Legendre 

Quadrature 

[1] 
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4.3 


The  Integration  of  the  Hull  Function  (P^.Q^) 


In  section  3.2,  were  defined  as: 


Qi 


L 

/ 

-L 


^(TK-l)^ 


sin 


(VKS)  dS  (4.27) 


cos 


and 


I;-; 

1-/N 

/‘  KUTK-nO 

l Q2J 

-L 

-T  C 

l J 

sin 


(VKO  dS 


(4.28) 


cos 


where  i = 1 is  for  pitch  and  i = 2 , for  heave.  Since  the  integrals 
in  the  braces  are  functions  of  K,  T and  S , they  can  be  easily  evaluated 
numerically  for  given  values  of  K,  x and  for  a given  section.  After  the 
sectional  integration,  integrations  along  the  length  of  the  hull  are  performed 
If  the  values  of  v,  T and  K are  given  and  if  we  let  the  inner  integrals 


in  (4.27)  and  (4.28)  be  f^(£)  , then: 


lpi' 

L 

sinaS 

= / MSH 

Qi 

-L  1 

cosaS 

v j 

MS 


1 - 1;  pitch 
i “ 2;  heave 


(4.29) 


where  u is  some  constant.  The  integral  of  the  form  (4.29)  can  be  evaluated 
using  the  concept  of  the  Filon-Trapezoidal  quadrature  [12].  The  basic 
idea  is  to  approximate  f.(S)  by  a linear  function  (a^S  + b^)  , say, 
between  the  intervals  and  S^+1  • Then  (P^.Q^)  can  be  approximated: 


K 

N C/+1 

sinaS 

| 1 

•=  l / (aX  + b)  < 

> 

K 

1=1  Si 

cosaS 
. - 

MS 


(4.30) 


wher 


e S,  - -L  and  SN+1 


V.  RESULTS  AND  DISCUSSIONS 

A computer  program  in  Fortran  IV,  based  on  the  analysis  of  this 
report,  has  been  developed,  and  calculations  have  been  performed  for  the 
sample  model  (Fig.  5)  for  several  Froude  numbers  using  the  IBM  360/370 
computer  at  the  MIT  Information  Processing  Center. 

As  input  data,  we  must  supply  the  offsets  of  each  section  of  the 
hull.  If  there  is  a parallel  middle  body,  we  need  to  give  only  the  offsets 
of  the  beginning  and  ending  sections  of  the  parallel  middle  body.  The 
offsets  of  the  remaining  sections  in  the  parallel  middle  body  can  be 
omitted.  For  Instance,  the  offsets  of  the  stations  9 and  11  are  sufficient 
to  take  care  of  the  parallel  middle  body  in  the  sample  model  (Fig.  5). 

The  most  time-consuming  computer  operations  occur  when  evaluating 
the  semi-infinite  Integrals,  especially  due  to  the  highly  oscillatory 
cosine  term  for  twin  hulls.  Numerical  convergence  of  these  integrals 
becomes  slower  as  the  Froude  number  increases. 

Calculations  were  made  for  Froude  numbers  0.0,  0.2,  and  0.4.  In 
Figures  6 through  11,  theoretical  heave  damping  coefficients  are  compared 
with  Lee's  experiments  [6].  In  Tables  3 through  5,  computer  outputs  of 
theoretical  pitch  damping  coefficients  are  listed.  It  is  to  be  noted  in 
particular  that  the  pitch  damping  coefficient  for  Froude  number  0.4  at 
low  frequencies  becomes  negative  (Table  5).  The  negative  damping  in  the 
present  study  is  not  easy  to  explain.  But  if  it  is  physically  realistic, 
it  implies  that  the  ship  will  be  unstable  in  pitch  at  high  speed.  The 
presence  of  negative  damping  was  noticed  for  an  oscillating  ellipsoid 
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near  a free  surface  by  Newman  [9]  and  by  Gerritsma,  Kerwin  and  Newman  [5] 
for  the  Series  60  hull  forms.  In  the  present  case,  where  the  hull  geometry 
is  such  that  zero  or  minimal  damping  can  occur  at  zero  speed,  one  antici- 
pates that  negative  damping  may  occur  sooner  than  was  observed  in  [5]  and  [9]. 

In  both  Figures  and  Tables,  it  is  to  be  noted  that  when  the  forward 
speed  effects  are  considered  there  are  peak  values  of  damping  coefficients 
at  critical  frequencies;  1.25  at  Fn  = 0.2  and  0.625  at  Fn  = 0.4 
for  which  T = 1/4  . Although  the  qualitative  agreement  is  good,  the 
theoretical  predictions  are  seen  to  be  somewhat  lower  in  the  higher 
frequency  range  and  higher  in  the  lower  frequency  range,  especially  near 
the  critical  frequencies  at  T = 1/4  . 

In  Figures  2 and  3,  thin-ship  and  modified  thin-ship  results  for 
two-dimensional  cylindrical  bodies  are  compared  with  Frank's  [4]  and 
Lee's  [6]  works.  Correction  factors  in  the  modified  thin-ship  approach 
may  vary  depending  on  the  hull  forms  and  the  frequencies.  At  high 
frequencies  the  modified  thin-ship  results  with  a correction  factor  of 
2 give  excessive  damping,  but  in  the  frequency  range  of  practical  importance 
this  modified  theory  agrees  better  with  Frank's  and  Lee's  results  than 
does  the  pure  thin-ship  theory. 

In  Figures  12  through  14,  damping  coefficients  are  plotted  against 
hull  distance  variations  for  Froude  numbers  0.0,  0.05  and  0.4  at  a 
fixed  non-dimensionalized  frequency  4.0.  In  order  to  investigate  the 
Brard  parameter  influence  on  damping  as  well,  Froude  numbers  are  chosen 
such  that  T = 0.0,  0.2  and  1.6  for  each  case.  In  Figure  12  it  is 
verified  that  heave  damping  coefficients  at  zero  forward  speed  oscillate 
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with  a period  of  A/2  and  the  mean  amplitude  of  oscillations  decays 
almost  like  (2vb)-1^2,  which  was  derived  in  (3.5)  by  the  asymptotic 
analysis.  As  the  hull  distance  and  Brard  parameter  increase,  the  oscil- 
latory phenomenon  dies  out  faster.  The  diminished  interference  effects 
at  higher  Brard  numbers  can  be  explained  physically  by  the  fact  that 
the  waves  in  this  case  are  confined  to  progressively  smaller  angles 
downstream.  Hence  the  waves  will  not  impinge  upon  the  opposite  hulls 
if  the  Brard  number  is  sufficiently  large.  The  minimum  hull  separation 
distance  at  which  the  interaction  effects  cease  can  be  predicted  by 
finding  the  generated  wave  angles  from  Fig.  1 in  [8]. 


15.36 

Midship  Section  (Station  10) 


Pig.  5 Dimensions  of  SWATH  Catamaran  for  the  Sample  Program 


STN 

X(I) 

EPS (I) 

TO  (I) 

R(I) 

1 

-261.00 

0.00 

41.68 

0.00 

2 

-250.00 

0.00 

41.68 

5.46 

3 

-240.60 

0.00 

41.68 

11.16 

4 

-230.00 

0.00 

41.68 

13.38 

5 

-220.00 

1.52 

41.68 

14.27 

6 

-209.00 

3.04 

41.68 

15.19 

7 

-198.00 

5.19 

41.68 

15.36 

8 

-176.00 

7.85 

41.68 

15.36 

9 

-154.00 

8.53 

41.68 

15.36 

10 

0.00 

8.53 

41.68 

15.36 

11 

154.00 

8.53 

41.68 

15.36 

12 

176.00 

9.72 

41.68 

15.29 

13 

198.00 

5.12 

41.68 

14.37 

14 

209.00 

2.87 

41.68 

13.48 

15 

220.00 

1.44 

41.68 

12.32 

16 

229.80 

0.00 

41.68 

10.62 

17 

239.50 

0.00 

41.68 

8.09 

18 

249.50 

0.00 

41.68 

4.61 

19 

259.00 

0.00 

41.68 

0.00 

Table  2.  Offaet  Data  for  the  Sample  Program 
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Froude  Number  « 0.0 


Table  3 


B55  = Pitch  Damping  Coefficient 
Damping 

" pVLVgt 

V = Displaced  Volume 

w * Encountering  Frequency  * u)  at  F =0.0 


weVL7g 

B5  5 

85  5 

Single  Hull 

Twin  Hulls 

0.5  ^ 

0.10851E-05 

0. 21689E-05 

1.0 

0.64874E-04 

0 . 12846E-03 

1.5 

0. 17057E-03 

0. 32428E-03 

2.0 

0. 10922E-03 

0 . 14321E-03 

2.5 

0 . 25516E-02 

0 . 15228E-02 

3.0 

0. 96172E-02 

0 . 37951E-02 

3.5 

0.17570E-01 

0. 13984E-C1 

4.0 

0. 20316E-01 

0.35164E-01 

4.5 

0. 18057E-01 

0 . 19299E-01 

5.0 

0 . 14401E-01 

0.33401E-02 

5.5 

0. 95581E-02 

0 . 14782E-01 

6.0 

0 . 54997E-02 

0 . 56590E-02 

6.5 

0. 28302E-02 

0 . 13068E-02 

Theoretical  Results  of  the  Pitch  Damping  Coefficients 
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Froude  Number  =0.2 

B55  = Pitch  Damping  Coefficient 

Damping 

” pvLVglr 

V = Displaced  Volume 

U)  = Encountering  Frequency 


uJUg 

B5  5 

Single  Hull 

Bs  5 

Twin  Hulls 

1.00 

0. 78276E-02 

0.39147E-02 

1.25 

0. 18242E-01 

0 . 23052E-01 

1.50 

0. 91238E-02 

0. 58606E-02 

2.00 

0 . 80448E-02 

0.57695E-02 

2 = 50 

0. 93458E-02 

0 . 70857E-02 

3.00 

0. 10716E-01 

0.13841E-01 

3.50 

0.11544E-01 

0.14518E-01 

4.00 

0 . 11595E-01 

0.71407E-02 

4.50 

0 . 10893E-01 

0. 12743E-01 

5.00 

0.97016E-02 

0. 13490E-01 

Table  4.  Theoretical  Results  of  the  Pitch  Damping  Coefficients 
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Bss  = Pitch  Damping  Coefficient 
Damp ing 

‘ pVLVir 

V = Displaced  Volume 

U)  = Encountering  Frequency 


CJ  \/l/ g‘ 

Bss 

Single  Hull 

Bs  5 

Twin  Hulls 

0.250 

-0.81562E-02 

0.66851E-03 

0.500 

-0. 64218E-02 

-0. 30374E-02 

0.625 

-0.12086E-02 

0 . 29785S-02 

0.750 

-0.24313E-02 

-0. 38668E-02 

1.000 

0. 20595E-02 

-0.21786E-02 

1.500 

0. 78046E-02 

0. 51413E-02 

2.000 

0. 84687E-02 

0.10173E-01 

2.500 

0.88749E-02 

0 . 11627E-01 

3.000 

0.88978E-02 

0.95881E-02 

3.500 

0.73347E-02 

0. 609C2E-02 

4.000 

0. 58270E-02 

0.46028E-02 

4.500 

0. 54560E-02 

0 . 54672E-02 

5.000 

0. 54788E-02 

0. 63667E-02 

5.500 

0. 52273E-02 

0.60154E-02 

6.000 

0.47029E-02 

0.49008E-02 

6.500 

0.39953E-02 

0 . 37948E-02 

7.000 

0. 31860E-02 

0.29722E-02 

7.500 

0. 24167E-02 

0. 23651E-02 

Table  5.  Theoretical  Results  of  the  Pitch  Damping  Coefficients 
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Heave  Damping  Coefficients  of  Single  Hull  MODCAT 


Heave  Damping  Coefficients  of  Twin  Hull  MODCAT  at 


□ 0.38 


Heave  Damping  Coefficients  of  Single  Hull  MODCAT 


Heave  Damping  Coefficients  of  Single  Hull  MODCAT 


1.8 


Fig.  11  Heave  Damping  Coefficients  of  Twin  Hull  MODCAT  at  F = 0.40 


Effect  of  Hull  Distance  Valrations  on  Damping 


Effect  of  Hull  Distance  Variations  on  Damping 
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APPENDIX  A 


INPUT  AND  OUTPUT 


INPUT  DATA 

DATA  CARD  I - FORMAT  (IIP,  2F10.2) 

NST  = The  number  of  offset  stations  of  the  hull. 

AL  = Length  at  the  waterline. 

B = Half  hull  separation  distance. 

DATA  CARDS  2 - FORMAT  (4F10.4) 

Data  cards  2 consist  of  NST  number  of  cards,  on  each  of 
which  the  values  of  X,  EPS,  TO,  R should  be  punched  according 
to  the  above  format. 

X(I)  = An  array  of  the  x-coordinatea  of  the  stations. 

EPS(I)  - An  array  of  the  half-beam  at  the  waterline. 

TO(I)  - An  array  of  the  distance  Tq  . (See  Fig.  5) 

R(I)  - An  array  of  the  radius  of  each  section. 

DATA  CARD  3 - FORMAT  (2110) 

NFR  “ The  number  of  Froude  numbers  to  be  tested. 

NOM  * The  number  of  non-dimensionalized  circular  frequencies 
to  be  tented. 

DATA  CARD(S)  4 - FORMAT  (8F10.4) 

FR(I)  » An  array  of  NFR  number  of  Froude  numbers. 

OM (I)  - An  array  of  NOM  number  of  the  non-dimensionalized 

circular  frequencies. 


L. 


- 49 


Besides  the  input  data  cards,  we  must  supply  the  roots  of 
the  Legendre  polynomials  and  the  weight  factors  for  the  Gauss- 
Legendre  quadrature  in  order  to  evaluate  the  semi-inf lnite 
integrals,  I^,  1^,  and  1^.  In  the  sample  program,  the  ten- 
point  Gauss-Legendre  quadrature  formula  io  used. 

OUTPUT 

Major  outputs  are  pitch  and  heave  damping  coefficients  of 
a single  and  twin  hulls.  However,  we  can  get  various  intermediate 
results  for  checking  purposes. 

DPI  (I)  = An  array  of  the  pitch  damping  coefficients  of  a single  hull. 

DP2(I)  = An  array  of  the  pitch  damping  coefficients  of  twin  hulls. 

DH1(I)  = An  array  of  the  heave  damping  coefficients  of  a single  hull. 

DH2(I)  = An  array  of  the  heave  damping  coefficients  of  twin  hulls. 
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PROGRAM  DESCRIPTIONS 

Main  Program 

The  main  program  handles  input  and  output,  calculates  various 
parameters  and  performs  six  integrals  in  (3.8),  (3.9),  and  (3.10) 
calling  subroutines  PIQI  and  ROOT.  Gauss-Chebyshev  quadrature  formula 
is  used  for  Ij  and  I3.  For  the  rest  of  the  integrals,  the  subroutine 
ROOT  is  called  in  order  to  determine  the  size  A for  the  sub-integrals. 
(See  Section  4.2.) 

Subroutine  PIQI 

This  subroutine  integrates  the  hull  function  in  (4.27) 

and  (4.28)  for  a given  section  first  and  then  integrates  along  the 
length  of  the  hull.  It  gives  the  values  of  (P_^2(K)  + Q^2(K)}  for 
given  values  of  X,  v,  and  T. 

Subroutine  ZETIN 

For  given  values  of  v,  T,  K and  for  a given  section  (i.e.,  e, 

T0,  r),  this  subroutine  evaluates  the  integral, 

-H  2 

yv(iK-l) 2 / h(£,OeV?(TK_1)  dC  (A.  1) 

-T 

by  using  Simpson's  rule  and  returns  the  resulting  value  through  the 
variable  Z.  The  y in  (A.l)  is  the  correction  factor  for  the  cylindrical 
part  explained  in  Section  2.3.  This  subroutine  is  used  for  both  (Pj.Qj) 
and  (P2,Q2). 
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Subroutine  ZETINS 


This  subroutine  evaluates  the  integral: 

/ h(£,C)CeVC(TK"1)2dC  (A. 2) 

-T 

by  using  Simpson's  rule,  which  is  similar  to  that  in  ZETIN.  The 
resulting  values  are  returned  through  the  variable  ZS  which  is  used 
only  for  (Pj.Q,). 

Subroutine  ROOT 

Thi.  subroutine  evaluates  the  zeros  of  cos  (2vb/(TK-l)4 - K2)  = 0 , 
from  which  we  get  a polynomial  equation: 

CtK-1)4  - K2  = (-^—)2  , (n  = 1,2,3 ) 

(A. 3) 

The  roots  of  this  polynomial  are  found  by  calling  the  IBM  Scientific 
Subroutine  Package  POLRT.  The  subroutine  ROOT  is  used  in  evaluating 
the  semi-infinite  integrals. 

Function  OS 

This  function  calculates  the  submerged  area  of  a given  section.  It 
is  used  in  evaluating  the  displacement  of  the  hull  in  the  main  program. 

Funf  tlon  F 

This  function  evaluates  the  integrand  in  (A.l)  for  the  purpose 
of  Simpson's  rule. 

Function  S 

This  function  evaluates  the  integrand  in  (A. 2)  for  the  purpose 


of  Simpson's  rule. 
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APPENDIX  C 
PROGRAM  LISTING 


1 


no  o non  non  o ooo 


T '•  r S PROG  Ml  ror^'IT’S  TF?  p ITCH  AND  PPAVF  DAMPING  COEPFTCTFNTS 
0 F A S I N G 1 r a s D TWIN  HIJtTS  SWATH  CC NP IG M n ATIO N 5 INCLUDING  PORWARD 
9 I = r 6 p»FpCTS,  8A6PD  ON  IFUIAN'S  THIN-SHIP  THfORY. 

Cln-NSTCN  TA'I  (10,  10)  ,040  (10,30) 

DTT1ISTC1  DM  1 ( 30)  , D'!2  ( K)  .DPI  ( 30)  .DP  2 (30) 

01“ =4  SIGN  7,  ( 10)  , ’4®TT  (10) 

0T1>  1STCN  PR  ( 11)  ,01  ( 10) 

CC  P ION  /GFOB/  NST,  NST,  PIS  (4  0)  , 1 (40)  , TO  (4  0)  , X (40) 

0C1BO1  /ONE/  FOOTS  (4)  ,FCCTT  (4) 

DATA  L/8/,  ,1/6/ 

SCOTS  OF  THE  10-POINT  (N  = 9)  LPGENDRF  POLYNOMIALS 
7(1)=. 14SH74 
7 (?)  =.  4 13  305 
7 (3)  -.870400 
7 (4)  =.  Pf  50*  3 
7 (8)  =.973006 


WFIGHT  factors  POO  thf  10-POINT  GAHSS-I.EGENDRK  QOADPATDRP 

WPTT  ( 1)  =. 206624 


W p T T ( ?)  = 
w ’ : ( D = 

.p^(4)  , 

WFTT  (9)  =■ 


260266 
2 1 00H6 
140461 
7666  7 1 


DO  100  1=1,6 
,1=T«8 
7 (.1)  =-->  (T) 


WCTT  ( .^ ) r W c T T (T) 

CC;TT1IUr 

tnph''  p»t.* 

0 1 0 (1,11) 

VST,  At.,  P 

V S ’ T (1,11) 

NO  r , A 1.  , 0 

->  ■ - C , 1.  1 

< (M  . -7  ; <T 

....  , - , , , 

(;."(’).• 

( n , 

p(  ' o (I.'-*) 

6P7, SOI 

wiT-r,  i8) 

■'  1,  M il 

nr  * 6 (1,16) 

(77(1)  ,'=1, 

•IPT'-rfB.  17) 

(=■  > (T)  ,T  = 1, 

f f>) 

R F A 7 (!.,1°) 

n* d) , t = i , 

* r-  T T p ( M , 1 0 ) 

(O'(T)  , T = 1 , 

>0») 

'fr  ( r 1 7 

,7713.  7) 

(D  , (r) 


11  =r  r “.\t  ( 17,  ' p A I'A: 

1 e V , * 1,  O . n = 

? ■ ! ■ "G  I I a T 


1 

“’A,*  ’I.’ 

•'ll!, T 

1 2 

nr)r»  v ^ I 

(4  = 10.4) 

1 1 

?r  ’■  v * r 

( 8 Y , 1 5 " N • , 

7Y  $ • 

i 

* J • //  ( 

'7.1 

ia 

°o  • *•  r 

(5T10) 

1 6 

eo‘,*r * 

(//6/,»|0. 

O9 

1 

6 v , • in. 

0* 

16 

»C I T 

( 3 9 1 ,1 . 4 ) 

17 

P ('  f W • *■ 

(•’‘•10,11) 

1 

v(. " M M 

(“V  1 1.4) 

10 

PCir  •• 

(P  FI  0.4) 

NS  T*  NST 

- 1 

SO^C  “OOPL  *OOC*.T  62?6'// 

MTTO'IS  = 1 , T 4/6 Y, 

IIP  VI.  =•  , '7.  2,  • (P^-'T)  •/ 

51  ICING  = • , ?7. 2 , ’ rFI'T)  ’////) 

( !)  * , 7v,  iTncji  ,7r,  * TO ' , 7X, 


’.rv/) 


r T = 1 . 1 4 1 o a ? 

r.  - 7 ? . I? 

T - 1 ^ 0 
; i-60 
.'2=10 


- —iOCIi 

•f  fuiJv 


",  ^ 4 

'■  »Ftad«c8c* 


MAIN0001 
BAIN0002 
BAIN0003 
B A IN0004 
BAIN0006 
HAIN0006 
BAIN 0007 
BAIN0008 
BAIN000O 
MAIN0010 
BAIN001 1 
BAIN0012 
BAIN001 3 
BATN0014 
BAIN 0016 
BA  I NOO  If, 
BAIN0017 
BA TNOO  18 
HA  I NOO 1 0 
BALN0020 
BAIN0021 
HAT  N0022 
HAIN0023 
MAIN0024 
BA  I N0028 
B A T N 00 26 
HAIN0027 
HAI N0028 

HAIN002O 
BA  I NOO  30 
HAIN0031 
BA  I NOO  32 
BAIN 0013 
HAIN0014 
BA  I NOO  36 
BAI NO 016 
BAIN 0037 
BATNOO 14 
i a r i i 

• .1  I • 0 4 ' 

1 A T 1 2 4 1 
BA  I N 0 0 4 2 
1 A I N 0 0 4 1 
BA  I N 0 0 4 4 
BAINOOuO 
B.ATN0046 
1AIN0047 
HAIN004R 
BAT  N 0040 
BAIN  0060 
BAT  NO 08 1 
BAI N 006  2 
BA  I 1008  1 
HAT N 0084 
BAT  NO 066 
B A I N 0086 
BA INO087 
BAIN0088 
BMN0060 
B A 1 11006 0 
BA  T N006  1 
BAI 10062 
BATNOOhl 
BATN0064 
B'  I N 0 0 6 6 
BAT  110066 
BA  T 40067 
BATN006R 
BAIN  1060 
BATNOO 70 
B A T N 0 0 7 1 
BATNOO 72 
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u u 


I 


UT 1 - n T /PI CAT (101) 
«i?=nr/?TC)\T  ( 10? ) 


ir>o 


?o 


« <■/ tCUUTTON  d'  TM*  CP1I-HULL  C I HP!.  T C "H  'NT  * 

> IN  CUBIC  PFFT  “"Cl  THC  HIV®*  DATA  • 

»•****•*♦**•*#***•*•**••••*••*••••*•*•••*•••••♦•• 

DISP=0. 

no  i5n  T=i,Nsr 
J - T ♦ 1 

F?S1=PPS  (T) 

PPS2=rPS  (J) 

R 1 = 8 ( I ) 

P 2 = ° (.1) 

TC  1 -TO  (D 
1 C 2 -TO  (,1) 

r i = x (i) 
x?*x  ( i) 

I)Et.T  = Y?-X  1 

nt1.V=SS  (”PS1  ,T01  , ?1)  *SS  (F.PS2,T02,R2) 

0 T S P " rI  SP  + o*’!. V*HF!. Y/2. 
whttf {",?'!)  ri;° 

POP  * t T (//6Y, •OISPLACENENT  =•.  F10.2,  • (CUPTC  FPPT)'//) 

»j-PTM  E S T O - T T T T NO  FACTORS  POP  PITCH  6 HPA.  VF  (S  I NOL  P HULL) 


HAIH0073 
HAIN0074 
HAIN0075 
H A IN0076 
HAIN0077 
HAIN0078 
HA  I N0079 
HAIN0080 
HA  I N 008 1 
HAIN0082 
HATN0083 
HAIN0084 
HAIN0085 
HAIN0086 
HAIN0087 
HAIH0088 
HAIN0089 
H1IN0090 
HA  I 9009 1 
HAIN0092 
HAIN0093 
HAIN0094 
HAIN009S 
HAIN0096 
HA  T N0097 
HAIN0098 
HAIN0099 
HAIH0100 
HAIN0101 
HAIN0102 
HAIN0103 
HA  X NO  104 
HAIN0105 
HAIN0106 
HAIN0107 
HA  I NO  108 
HAIN0109 
HAIN0110 
NAIN0111 
r.  AIN01  12 
HA  T NO  1 1 3 
H A T NO  1 14 
HA  I NO  1 1 5 
HAIN01 16 
HAIN0117 
H A I NO  118 
HAIN01 19 
HAINO 120 
HAIN0121 
HAIN0122 
HAIN0123 
HAIN0124 
NAIN0126 
HAIN0126 
HA  I NO  1 27 
HAIN012B 
HA  I NO  129 
HAIN01 30 
HA IN0131 
HAIN0132 
HAIN0133 
HAIN0134 
HAIN01 35 
HAIN0136 
HAIN0137 
HATN0138 
HATNO 1 39 
H A T N 0 1 4 0 
HAt  NO  IN  1 
HATN0142 
HA  I NO  1 4 1 
HATH  I)  144 


no  on 


IP  (Z1.IT.0.)  11=0. 

A FG  = 2. *SNU*?*S0RT(Z1) 

C 

C7  I L PIC  I (XKT,A  Nil,  oTA0,  F rc,  >1P0) 

c 

C I NT  p 3HE  T T ATS  RESULTS  CAN  PE  PRINTED  IP  DFSIRED 

C HRIT1(P,23)  TKI,*OG,Pi>C«fiPQ 
C 21  PORHAT  (2F10.4,  1I.2M6.5) 

C 

FPO=-PPQ 

PPO=-MPC 

C 

PI  1 = PI  1 4PP0 

PI2  = PI?«rFQ*2.*  (COS  (ARG) ♦ 1. ) 

HI1=HI1*HF0 

RI2=“T2+HFp*2.* (COS (ARG) ♦ 1. ) 

1400  CCNTINUP 

c 

DF1(.1)=((PI1*WI1)*  (-2.  «RN*A  NII/PT)  )/PNON 
DH  1 (.1)  = ( (PI  1*HT  1)  * (-2.  « VI*I  NII/PT)  ) /HtlON 
11P2  (.1)  = ( (PT2*WT  1 ) * (-WN«AN'I/PI)  )/PN0N 
CH2  (J) = ( (PT2*WT 1)  • (-RN*f NU/PI) )/HNON 

GO  TO  4 10 

ISO  ^ C V T I N 0 F 

RK 1 , RK2, PK3, RPR : ROOTS  CP  (TAU*K- 1 . ) * • 4 -K ** 2=0. 

RK 1=0. 

R F 2 = 0. 

RK  1 - 0 . 

PKU=1. 

T A 1 = Sop  1 (1.-4. • UTAH) 

T r ? = SCPR  ( 1.  ♦<!.  *:'T'  0) 

Tf  1=2.»CTM)-1. 

T"2  = 2.»iT  M! ♦ 1. 

T Al  » 2.  *r  T A 

r * i--  ( ; •;-■•.'■)/ : \ j 
FK4  = (1:  2*11  2)  /T.**" 

I f (IT,-  U.  10.  3)  GC  TO  5CC 

TA  1--S0PT  (1.-4.  »PTAU) 

IP  1=  (1  1-11  1)  /TAP 

R«2=  (IF  1<TA1)/TAD 
S00  CCNTINO' 

WRITE  (V,2H  I<K1,1K2,RK 
21  FC"M'1  (//r>X,'K1=»,  PH.f.«X,,K2*',E13.b, 

1 4X  , 'N  3='  , ?.  13.  f , 4X,  *K4  = • , E 11.fi/) 

C 

P H2'0. 
f 212=0. 

H1I?=C. 

1121  2 = 0. 

C 

II  (ITA  U.FO.  1)  GO  TO  S5C 

C 

c 

c 

C » G < T A 0 < 0 . 2 S : PINTTf  I N T P(i  1 A L PPCE  K2  TC  Kj  • 

C • USING  '(,2-IOINT  GAUSS-CURT  Y3IIE  V QUADRATURE  * 

C ••*•*•*«•«•*»<•*•••**•••«*••*•**»••*•»***••**••*• 

C 

A 1=  (3K3-RP2)  /2. 

E 1 = (RK3«FF2)/2. 

C 

DC  6)0  J2=1,»02 

XET=COS (?tr • T (2»J2- 1) *F I/FLOAT  (2*H02) ) 
AKI=A1*XKI»P1 


HAIN0145 
HA  I NO  1 46 
II A I NO  I 47 
HAIN014R 
HAIN0149 
HA  I NO  ISO 
HAIN01S1 
HAIN0152 
HAIN01S3 
HAIN01S4 
HATN01S5 
HAIN01S6 
HA  I NO  I 67 
HAIN01SB 
HA  I NO  1 S9 
HAINO 160 
HAIN0161 
HAINO 162 
HAIN0163 
HAINO  164 
HAIN016S 
HAIN0166 
HAINO 167 
HA  I NO  160 
HAIN0169 
HAIR0170 
HAIN0171 
HAIN0H2 
HA  I NO  1 7 3 
HAIN0174 
HAIN0175 
HAIN0176 
HAIN0177 
H A I NO  1 78 
HATN0179 
HAIN0100 
HAIN0181 
H ; I H 0 1 P 2 
* A I j;  1 H 1 
M A 1 1.  <11  8 4 
HAIN018S 
H A I N 0 1 86 
HA  I NO  1(17 
HAIN01H0 
HAIN0189 
H A I NO  1 90 
HAIN0191 
HAIK0192 
HAIN019J 
HA  I NO  194 
HAIN019S 
HAI  NO  19»i 
HAIN0197 
HAI N 0198 
HAIN0190 
H.AIN020  ) 
HA IN020  1 
HAIN0202 
H A I NO  20  3 
HAIN0204 
HAIN020S 
HAIN02O6 
HA  I N 020 7 
H A I NO  20  3 
HAIN0209 
HAI N02 10 
HAT  N02 1 1 
HAIN0212 
HAIN0213 
HA  I N0214 
HAIN021S 
HAIN0216 
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on  nonoonnnoon 


L 


C 

Z 1= PTAO**2*S0RT  ( (RK4-AKJ)  * (AKI-RK1)  ) 

Z2  = (BTf  U*AKl-1.)  »*4 
23  = 7, 2-AKI**2 
If  (ZJ.1T.Q.)  73=0. 

Z4  = C0S  (2.  • ANl'»3*S3RT  (Z3)  ) 

C 

SGN=  1. 

SCNX=ETAC*AKI-1. 

If  (SGNX.IT.O.)  SC.N  = -1. 

C 

27  1=Z2*SGN/7  1 
ZZ2=ZZ1«  (2. *74  + 2.) 

C 

CALL  Ficr  (AKI, AND, BTAC, PFC, HPC) 

C 

FI  1 = F I 1*PFC*7.7.1 
FI2=FI2*EF5*7Z2 
HI  1=HI1*HI0*ZZ1 
HI2=HI2«HF0*ZZ2 
600  CCNTINLE 
C 

F1I2=PI1*VI2 

F2I 2=  FI  2 *MI 2 
H 1 1 2 = HI  1 *fc  r 2 

h;i  ;=Hi2»wi2 

c 

C ««♦««.♦«♦«•••»»*«*«*•••*«•*********•••«•*•••♦••••• 

C • II  A U=  2 : INI FG RATION  FROM  *-•  INPINITY  TO  K1  • 

C • IT  AU=  3 : INT:Gr'ATICN  F FOB  INPINITY  TO  A3  * 

C 

C 

560  CCSTINUf 
S 6 B P 1 = 0 . 

su;ir2=c. 

El  All  1 = 0. 

SH*M2=0. 

» - n 

X < > • I “•  1 

ii  ('.  ; . i ■ . 1 ) x v ( = :i  k 1 

ro  650  ITFR=1,Nr,7 

c 

CALL  ROCI  (ITF3, 37Ar/,ANU,P) 

C 

IF  ( 1,nOT  R (1)  .Lr.XXX.ANC.RCOTI  (1)  .50.0.)  CXX  = FCOTP(I) 

IF  (l-CC  13  (2)  -Lr.X'X.  A NT.  FOOT!  (2)  . FJ.O.)  CXY=RCOTP(2) 

IF  (ROGIB  (3)  .LT.  < X X . A HD . ROOT  I (3) .EC.O.)  CXX  = PCCTP(3) 

IF  (RCCTF  (4)  .IT. XXX. AFC. FOOT I (4)  «E0*0.)  CXX=ROOTR(4) 

E I = 50rT  (VXX-CXX) 

XI  !:  = "!-. If 

1ST "RTF  0TITc  RESULTS  CAN  PE  PR  T NT  E C IP  OESIRFC 

5 F T T " (',30)  1 TCR 

30  ter.-  r ( ix , « it? ;<  = • , 13/} 

Afir  (B,1)  (ROO.-R  ( I 1)  , 11=  1,4)  , (RCCTI  (I  1)  , I 1 = 1,4) 

31  FCPBfT  (MX,  <i.  1 = • ,X13.F  . * X,  • R2=*.  F:13.6,  6X  , • PJ  = « , P13.6.6X,  • R4»*  , 

1 PI  1.6/HX,  M 1=*  , 111  3. 6,  *.  X,'  1 2 = *,  El  J.  6,6X  , • 13**,  E 1.1.6,  5X,  • 14**, 

2 '•13.6/) 

WF1  I P («  , 3 2)  A A , "B,a  X 

3  2 FCl'BAl  (RX,  • AA  = * , El  1.6, 4X,  • FP*  • , E 13.6.4X,  *CXX=F  ,R13.  6/) 

PI  1 = 0. 

FI  2 = 0. 

H 11  = 0. 

>112  = 0. 

to-iciNT  r.\,r ■•;-[..=>;pmppf  coa  r./A  tore  pen  a oiv'n  intprval 

CO  7)0  .11=1,10 
x a ',*!»=(  x r f • z t.T  i )♦•■•?♦:!)/  i . 

» I VEIT  (.13) 


F1AIN0217 

(1AIN021  8 

BAIN0219 

B AI N 0220 

BA.IN0221 

BAIN0222 

BAI N0223 

BAIN0224 

HAIN0225 

BAIN0226 

BAIN0227 

BA  IN0228 

BAIN 0229 

BAIN0230 

BAIN0231 

BAIN0232 

BAIN0213 

BAIN02  34 

BArN0235 

BAIN0236 

BAIN02J7 

BAIN0238 

BAIN0239 

BAIN0240 

BAIH0241 

BA.IN0242 

BAIN0243 

BAIN0244 

NAIN024S 

BAIN0246 

HAIN0247 

BAIN0248 

BA  I NO  2 4 9 

BAIN02S0 

BAI N0251 

HAIN0252 

BAIN0251 

BAIN0234 

BAIN0255 

BA.  T N 0 256 

B A I N'  0 2 3 7 

BA IN1258 

KAIN02S9 

BAIN0260 

BAI N026 1 

BAIN0262- 

HAIN0263 

BAI N0264 

BA  I NO  265 

BATN0266 

BA IN0267 

BAIN026H 

BAIN0269 

BAI NO  270 

BAIN027 1 

BAIN0272 

FIAIM027J 

BAIN0274 

BAIN0275 

BA  I N0276 

BAIN027  J 

BAI N027M 

BAIN0279 

BAI N02H0 

BAI N02H1 

BAI N0282 

BAI N071  I 

MAIN0784 

BAI N02OS 

B A I N 0206 

BAI  NO  >->7 

BAI N07O8 


1 

; 


- able  to  DUO  does  ool 
lull-/  It  Mo  xeptotluctioo 
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n o 


■ ' IS  BBS  BiiES  NOT 

r'trlfwm  maam 


* ! t,  LEGIBLE 

A*I=PK1-IPBAR»*2 

If  (ITAU.EQ.3)  AKI=RK3-XKPAB»»2 

IP  (ITAU.F9.3)  GC  TO  71C 

Z 1-f TAU»*2»SCRT ( (3KU-AK I) » (BK  J-AKI) * (RX2-AKI) ) 

GO  TO  72 C 

710  Z1  = rTA'l*SCRT(((nT10*AKI-1.)**2«AKI)*(RIC4-AKI)) 

720  Z2=  <rTA’J»AKI- 1.)  «*4 

£GN= 1. 

SGNX=ETA0*AXI-1. 

If  (SGNX.LT.O.)  SGN=-1. 

Z3=Z2-APT**2 

23  SHOUID  MOT  RE  NEGATIVE  VALUE. 

HOWEVER,  NUMERICAL  EPPCE  MAY  GIVE  VERY  SHALL  •-*  VALUE. 

IF  (7.3.LT.0.)  7 3=0. 

Z4=COS  (2.  »ANi:»3«SQRT  (2  3)  ) 

ZZ  1=72* SG  N/Z 1 
ZZ2=ZZ 1 ♦ (2.*ZU*2.) 

CALL  PIQI  (AKI,ANU,ETAU.PPQ,HPC) 

PI1=PI1«XIM*WI*PP0*ZZ1 
FI2=PI2*Xlw*yi»FF3»ZZ2 
RI1=HI1 ♦XIN«WI*HFg«7Z 1 
HI2-HI2*XIN*WI*HFJ*ZZ2 
700  CONTINUE 

SUMP1=SUFF1*FI1 

SLMP2=SUMF2*PT? 

SUMH1-SUPK1*HI1 
SUM!I2  = SU«H2*!II2 

EEOPP=AES (PI  1/SUM  P 1 ) 

ER0P'1=AES  (HI  1/SUEH1) 

RETT"  («,13)  PI1,S'I'*ol,EP0RP,PT2,SU*r2 
rsr  * ‘ C I ' T t = T 5 C.'f  IE  1 'I-.’TfC  TP  DrSIT-rr 

33  ec. >• ■-  •:  i-  x,  • pi  i-  •,  ••  i ?. h,  n,  • su-t  i = * ,i  1).'-,4<,*eropp=',  = io.3, 

1 4X,  •ri2-*,  - 1 ’.3,4*,  •ZU-p.’-'.E  13.'./) 

WRIT-  {'•,34)  HI  I.SHIM.E  HORH,'II2,SUPH2 

34  FQRM.'T  (EX,  • HI  1 = •,  El  3.6, 4 X,  • SUMH  ! = • , E 13 . 6, 4 X , • F ROHH=  • , E 1 0.  3 , 

1 4X , • H I 2 = 1 , E 1 3.6, 4X, •SUMH2=* ,E 13. 6/) 


A A = P n 

rf  (EROFF.IT.  10.  F-05.AU  .EPORH.LT.  10.E-05)  GC  TO  800 
630  CONTINUE 
800  CONTINUE 

PIII-SUBPI 
P 2 I 1 =SU  HP  2 
H1T1  = SMPM 
H2I 1=SUMH2 


TMDO.  : INTEGRATION  PPOM  K4  TO  INPINITY 


S0MP1=0. 

SUMP2=C. 

SUM  >11=0. 

SUM H 2=0. 

A A*0. 

DC  8r»0  ITER*!, NT, 2 
C A 1 1 PCCT  (ITER,r'TAU,3RU,P) 

IP  pCOTE  (1)  .or. RN4. A fC.NCOTr ( 1)  . EC. 0.)  CXX  = ROOTR(1) 
If  ( ROOT  E (2)  .GI.PK4. AND. ROOT I (2) . EC»  C . ) CXX  = FOOTR(2) 
IE  (ROOT  R (3)  .Gr.FT4.ANC.RC0H  (3)  . EC  • 0.  ) CXX  = FCOTF(3) 
If  ( ROOT  f (4)  .GT.RX4.ANr.R00TI  (4)  . EC  . C.  ) C X X = ACOT  R (4) 


BAIN0289 
HAIN0290 
HAIN0291 
HAIN0292 
HAIN0293 
HAIN0294 
H A I N0295 
HAIN0296 
HAIN0297 
HAIN0298 
MAI N0299 
HAIN0300 
HA  I NO  301 
MAIN0302 
MAI N0303 
MAIN0104 
HAINO  103 
MAIN0306 
HAIN0307 
HA  I NO  308 
HAINO  309 
MAINOJ10 
HAINO 31 1 
HA  I NO  11 2 
HAINO J1 J 
HA  I NO  3 1 4 
HAIN031S 
HAINOJ IN 
HA  I NO  3 1 7 
HAIN0118 
BA  I NO  J 19 
HAIN0J20 
HA  I NO  32  1 
HAINO 122 
HAIN032 3 
HAINO  324 
HAT  NO  323 

r,t. tno  i.*» 
H At N 0 3 7 7 
HA  I NO  I M 
HAINO  329 
HAIN0110 
H A 1 11 0 3 3 1 
PAIN0332 
HAIN0331 
HA  I NO  3 34 
HAIN0335 
HAIN0336 
HAIN0337 
HATN0338 
HAIN0339 
HAINO  140 
HAINO J41 
HAIN0342 
HAIN0343 
HAINO 344 
HAINO 143 
MAI  NO  343 
HAINO 34 / 
HAIN034H 
H A I N 0 3 4 9 
MAI  NO  130 
HAINO  IS1 
HAINO  132 
HAI  NO  13  3 
HAINO  194 
HAINO  133 
MAIN019h 
HAINO  137 
HAIN033H 
HAIN0339 
HAIN 0380 


u u 


^ IU.IMIB.pl. 


m AVAILABLE  TO  DOG  GOES  KOI 

(OLE  PRODUCTION 


- - ^ *«t  r".i  i \e  1 r 

i I s t; 


LflJ 


c 

EE=SORT (CXX-FA4) 

xin=?e-aa 

c 

C INTF!RNPCI?TP  RPSOLTS  CAR  PE  PRINTED  IP  DESIRED 

C SPITZ  (E,UO)  ITER 

C 40  FCRRAT  (3X.  • ITER**, 13/) 

C SPITE  (S.41)  (ROOTR  (12)  ,12=1,4) , (RCCTI  (12) .12*1,4) 

C 41  PCRRAT  |8X , • R1  = • , 213. 6. i X. • R2=» , El 3.6.5X, • R3  = * ,E1 3.6.5X, • H4  = * , 
C 1 EU.6/8X,,I1=,,E13.6,?X,*I2=,,E1J.b.5x,,I3*,.E13.6,5X,,I4='. 

C 2 E13.6/) 

C SPITE  (4,42)  AA.3B.CXI 

C 42  PCRRAT  (8X,1  AA  = * . E13.6.4X,  • EE=»,M3-6,4X,  •CXX=*  ,E13.6/) 

C 

PI  1 = 0. 

FI  2=0 . 

R 1 1=  0. 

1112=0. 

....10-POINT  GAIISS-LEGFNCRF  CUADRITORE  PCR  A GIVEN  INTERVAL 
DC  90C  *14=1,10 
XPPAR=(XIN»7(J4)«AA*EP)/2. 

S I = WErT  (J4) 

ANI=RK4*XHBAR*#2 
IF  (ITAt.PO.  3)  GO  TO  9 1C 

Z 1=  27AU  • •2*SCR  r ( (Atr-RP  1)  • (AKI-RR2)  * (AKI-FR3)  ) 

GC  TO  420 

910  Z1=PTA0*SCpT  ( ( (0TAO=AEI-1.)  ••2*ARI) • (AKI-RK3)) 

92C  Z2=  (ITA0«AKI-1.) *»4 
C 

SGN=  1. 

SG  NT  = F T A U*  A K I- 1. 

IP  (SGNX.LT.O.)  SGN  = - 1 • 

C 

Z3=Z?-AKI**2 
IF  (23.LT.O.)  7,3  = 0. 

Z4=C0S <2.»  A NO* 9* SORT  (23) ) 

771  = 7. ’•S”1\/T1 
’".’=  ' ' 1 • (.  . *’4*  >.  ) 

c 

c.il  uci  (M,.'n,-'iiu,nc,ii?i) 

c 

Fii  = rn*xiN««i»rpo*7.zi 
FI2  = °I?«XIN»WI»PP3»7.22 
111  1 = II I 1 ♦ X I *1*  * I •HPO*’tZ  1 
HI?=II2*XIN*SI*HrC*772 
9CC  CONTINUE 

c 

SUNP1*S0PP1*PI1 
SU.1P2  = SCPr2*Fl2 
S "J*  '1 1 =51  E !l  1 ♦TIT  1 
SUNIl2=SDPH2*III2 
C 

E ROP  P=  A F S (411/30741) 

FPOfl'l=AF£  (HI  1/SUEH1) 

C 

C INTFFPFCIfTE  RESULTS  CAN  RE  PRINTED  IP  DESIRED 

C SPITE  (.*,50)  PI1  .S'lRFI,  F BCRP , PI2 , SOU P2 

C 50  FCR.V.I  (EX,  • FI  1='  . '1  3.6.  4X,  'SURP1**  , F 1 3.6.4X,'  ERORP  = » .E10.3, 

C 1 4X,fPT2  = ,,c*13.P,4X,,S  II  4P2=I ,313,6/) 

C SPITE  (R.51)  HI1,SII*1H1,IFORH,HI2,SURP2 

C 51  PC  RE  AT  (6  X , • II I 1 = * , F.  1 3. 6 , 4 X , • SO  E H 1 = * , 4 1 1.  b,  4 X , ' E RORH=  • , P 1 0.  3, 

C 1 4X,  *11 12=' ,413.1,41,  'SUN II 2**,  H 1.6/) 

C 

A A » F !• 

IE  (Ff'OPF.LT.  10.  »-05,  ANC.EPORH.lt.  IO.p-05)  GC  TC  950 

050  CONTINUE 
950  CfhTtNUP 

c 

[ IT  1=  SO  E P 1 
Pi  1 3-3UEP2 
II  It  3*511411 1 


- 60  - 


m 


HAIN0361 

HA  I NO  36  2 

H A I NO  36  3 

RA I NO  164 

H A I NO  365 

BAINO  366 

RAINO  3b7 

HAIN0368 

RAINO  369 

HAIN0370 

RAI N037 1 

HAIN0372 

RAIN0373 

RAIN0374 

RAI NO  375 

HAIN0376 

HA  I NO  377 

HAIN0379 

RAIN0379 

RAINO  380 

RAIN0381 

HAIN0382 

HAIN03H3 

RAINO  384 

HAIN0385 

HA I NO  386 

RAINO  387 

HAIN0388 

RAIN0389 

RAIN0390 

HAIN0391 

RAIN0392 

RAIN0393 

RAINO  394 

RAIN0395 

RAINO  396 

RAINO 397 

HI  I NO  399 

RATN0399 

HATN0410 

r.  A I NO  40  1 

RAI N040? 

RA I N040  3 

RAIN0404 

RAIN0405 

HAIN0406 

HAIN0407 

RAIN0408 

RAIN0409 

HAIN0410 

RATN041 1 

HAIN0412 

HAIN041 3 

P.AIN0414 

RAIN0415 

HAIN0416 

RAIN0417 

HAIN0418 

HAIN0419 

RAINO 4 20 

HAIN0421 

RAIN0472 

R A I N 0 4 2 I 

RAIN0424 

RA I N0425 

H A I N 04  26 

RF  t N 0 4 2 7 

RAIN04.MI 

P.AIN0479 

HA  t NO 4 10 

HA1N04  I 1 

RATN04  12 


aSI3=SUP>H2 

c 

DF1  (.1)  - ((MI1*P1I2*P1I3)*(-2.  *WN»ANn/PI)  ) /PNON 
0H1  ( J)  = ( (H1I  1*il1T2*H1I3)  * (-2.*liN*ANU/PI)  ) /HNC  N 
CF2(J)=(|F2I1*P?I2*P2I3)»  (- 9 N • A HU/F I ) ) /P NCN 
0H2  (J)  = | (M 21  1*II2I2*H2I3)  * (- N N* A NU/P I ) )/HNON 
C 

MO  CCNTINtJF 
250  CONTINUE 

SPITE  (“.52)  (OB  (J)  ,DF1  (J)  , EP2  (.1)  ,DH1  (J)  ,DH2  (J)  ,J=1,NOfl) 

52  FCEBAT  (///5X, •CNEGA'.  1CX.  • DM* , 1 3X  , • DP2  • , 1 3X  . • CH  1 • ,13X, *DH2 •// 
1 (51, P5. 3, IX, 4*16. 5)) 

200  CCNTINUF 
STOP 
FDD 


I 


BAIN0433 
BAIN0434 
BATN0435 
BAIN0436 
HAIN0437 
BAIN0438 
B AIN0439 
BAIN0440 
BAIN044 1 
BA  I N 04N2 
BAXN0443 
BAIN0444 
BAIN0445 
BAIN0446 
BAIN0447 


‘ 


SUB  PQUTI  N 5 PTOr(AK,flNC,1AU,PPO,HPO) 

PIQI0O01 

CC'1  /OH.'./  NSI,NSI,iFS(4C)  ,R(U0)  ,TC(4(1)  ,X  (40) 

PTOI0O02 

r i - o . 

P 1 01 000  3 

r : - o . 

PI0T0004 

C 1 = 0. 

PI  01 000  5 

C2=0. 

P I 01 0006 

c 

PIOIC007 

DC  1 1=1, AST 

piQioooa 

c. . . 

...NSI  IS  THF  NO.  OP  THF  INTEGRATION  I NTERV AL ( NSI=NST- 1) 

PIQI0009 

J«I«  1 

PIOI00  10 

FFS1=EPS (I) 

PIQI0011 

F FS 2=  E PS  (.1) 

PIQI0012 

R1  = P (I) 

PI0I0013 

R 7 = R (.1) 

PI  01 00  1 4 

1C  1 = TC ( I ) 

pioioo  is 

TC2  = TO  (.1) 

PI0I0016 

H H 1 = SC  PT  (M*«2-»P51**2) 

PIQT0017 

H H 2 = SC R T (P2»»2-:PS2**2) 

PIQIOOId 

H 1*T01-HM 1 

PIQI00 19 

H2=T02-HH2 

PI0I0020 

X 1 = X (I) 

P I QI 002  1 

X2  = X (.7) 

PI  01 0022 

DFLX=X2-X 1 

PI QI 0023 

tt  = PUU»  (T  AU*AK  - 1.)  »»2 

PIQI0024 

ttc=nc»rt 

PIOI0025 

I RUK  = PtlU*»K 

PIQI0026 

c 

PIQI0027 

11  1*!>NUK»X1 

PIQI002B 

XX2=0NIIK»X2 

PIQI0029 

c 

, 

PIQM030 

CALL  ZFTIN  (DD, PFS1, TCI, R 1,21) 

PIQIOOJ 1 

CAIL  7F.TIK  (CD,EPS2,TC2,R2,Z2) 

PI0IO032 

c 

PIOIOO  3 1 

call  inns  (rn,iPsi,Tci,Ri,Hi,zsi) 

P1QI0034 

CALL  ZtTINS  (3D, - PS2 . TO J , R2 , H2, 2S2) 

PI0IOO35 

c 

PIQI0036 

nnn 


i 


W H 1 = -H  1 * I C 
WW2=-H2»CD 
C 

IF  (WW1.1T.-S0.)  EP1=0. 

IF  (WW1. IT. -50.)  GO  TO  10 
FI 1 = CXP  ("hi) 

c 

10  IF  (4W2.IT.-50.)  EP2=C. 

IF  (4W2.II .-50.)  GO  TO  20 
FI 2=EXP  (WS2) 

20  CONTINUE 
C 

PH 1=PFS1*FP1-Z 1 
PH2=EP52»FP2-Z2 
C 

P F A 1-Z  5 1 

PF  A 2 = ZS2 

C 

AA1=((-4W1*1.)/DCD)«FF1-(1./CEC) 
AA2=((-W12*1.)  /ODD)  *FP2-  (1./DDC) 

P F B 1=*F  P S 1 • A A 1 
FFF.2=EPS2«/F2 
C 

FrC1-X1*FHl 

FPC2-X2*FF2 

C 

A I A I ■=  (PF A2-FPA1)  /DELX 
AFBI=(FFE2-FPl!l)/DELX 
A FCI  = (PFC2-FFO 1) /DELX 
C 

BFA I*FFA 1-APAI»X1 
BFni=FPE1-APII«X1 
FFCI-PPC  1-f  HCT"! 1 
C 

IF  (BNUN.FQ.O.)  GO  TC  SC 
C 

FFS1=(rC|XX?)-PC(XX1))/E'II!K 
F FC  1=  (F  S (XX2)  - "S  (YX1)  ) / * SUP. 
r 

F 3 1 - . * \ I v v ? ) - •>  1 S ( y < 1) 

FCC1  =Ci  S (XX  1)  -COS  ( X X 2 ) 

C 

FFS2= (PS (XX2) -PS (CXI) ) /ENUK**2 
PFC2=  (FC  (XX21-FC  (XX  1)  ) /E  NTJK**2 
C 

F £ S2=  (COS  (XX 1)-COS  ( X X 2)  )/BNUK 
F C C 2 = (SIN  (XX2) -SIN (XX 1) )/INUK 
C 

A H I-  (FC2-FH  1)  / CPI. X 
DPT=PH1-AHI*X1 
C 

FI 1 = - ( ( A F A I ♦ P TP  I) * PP  S 1 ♦ (FPAI  + 8PPI) *FS51) 
1 - (APCT*?rS2*3PCt,»FSE2) 

C c i = (?r.rT*;r,i)*F=,ci«(iPir«FPFi)*Pcci 
1 - (AFC  I«F?C2*:i?CI*?CC2) 

F F 2 * A H 1 »FFS2«r  1 1 • F S S 2 
CC2*AHr»FFC2*P*II  • PCC2 
C 

P1*P1«PF1 

c 1 -C  1*<>C  1 
I2=P2*PF2 
C 2»C2*C02 

C 

GC  TO  1 
SO  CONTINUF 

CAST  FOB  r’(UK  = 0. 

1 1«  r?=o. 

Cl  =01- (FPC 1 *PPC2) ‘DEI X/2. 
g 2 = Q2* ( FH 1 ♦ ®M2) • D8LX/ 2« 

C 

1 CONTINUE 


PIQI0017 

PIQI0038 

PI0I0039 

Pigroowo 

PIQI004 1 
PIQI0042 
PIQI0043 
PI QI 0044 
PIQI004S 
PIQI0046 
PIOI0047 
PI0I0048 
PIQI0049 
PIQI0050 
PIQI0051 
PIQI0052 
PIQI0053 
PI QI 0054 
PIQI00S5 
PIQI0Q56 
PIQI00S7 
PI QI0058 
PIQI00S9 
PI QI 0060 
PIQIOOB 1 
PIQT0062 
PI QI 00b  3 
PIQI0064 
PI QI 006  5 
PIQI0066 
PIQI0067 
PIQI0068 
PIQI0069 
PIQI0070 
PIOI0071 
PI QI 0072 
PI Q 10073 
PI  QI 0 07 4 
PIQT0C75 
ri  Qior.76 
PI QI 00  7 7 
P10T0078 
PIQI0079 
PI  010080 
PIQI0091 
PIOI0082 
PIQIOOH  3 
PI QI 0084 
PIQIOOBS 
PI QI 0096 
PI QI 0097 
PIQI0068 
PIQT0089 
PIQI0090 
PIQI0091 
PICI0092 
PIQI009 1 
PI  QI 0094 
PI QI 0045 
PI QI 009b 
PIQI0C97 
PIQI009H 
PI QI 0099 
PIOIOIOO 
PIQT0101 
PIQI010? 
PIOI 010  1 
PIQIO 104 
PIQI0105 
PIQI0106 
PIQI0107 
PI QI 0 1 0 8 
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SO  F BOUT  INF  Z FTI N ( DO, P PS , T , F , Z) 

ZETNOOO 1 

IMH.EQ.O.)  2 = 0. 

2ETN0002 

IP(R.EC.C-)  PETURN 

Z ETNOOO  3 

MiR  = 10 

ZETN0004 

EIV=2.*PICAT  (NON) 

ZETNOOOS 

A=  F.PS/P 

ZETN0006 

CANNA=2. 

ZETN0007 

A1PA=ASTK (A) 

ZETNOOOS 

E=F*COS (ALPA) 

ZETN0009 

TA=P-T 

ZETN0010 

TE=-T-R 

ZETN001 1 

F,T  AFP  POSITIVE  ( TRE  EACTOUS  ANE  DEPTH  TO  THE  AXIS) 

ZETN001 2 

H= (TA-TF) /DIV 

ZETN001 3 

I = TB 

ZETN001 4 

2 = 0. 

ZPTNOO 15 

EC  1 i=i,Nnn 

ZETN0016 

X H = X ♦ H 

ZETN0017 

XHH=  X*2. »H 

Z ETNOO 1 8 

2=7*  !1/3.  » (F  (DD.T.R.X)  *4.  *F(CE,T,R,XB)  *F (DC.T.R.XHH)  ) 

ZFTN0019 

X=X*2.*H 

ZETN002C 

z=oanr; *2»nD 

ZETN002 1 

RETURN 

ZFTN0022 

EAC 

ZETN0023 

; 


S I)  r ROOT  I N E ZFTINS  (DC  , E I S ,T0,  R , H . ZS) 

ZTNSOOO 1 

C 

tt  = N\l*  (TAI)»P-  1 ) * * 2 

ZTNS0002 

IE  (;'.rC.U.)  2 S = 0 . 

2TNS0003 

IF  (R.iC-C.)  RETURN 

ZTNS0004 

N=  10 

ZTNS0005 

EIV=2.*EICAT  (N) 

2TNS0006 

T = TO*  F 

2TNS0007 

F = (T-HJ/EIV 

ZTNS0C08 

C..  . 

..•M  IS  THE  SFACINJ  Of  TEE  INTEGRATION 

ZTNS0009 

X = -T 

2TNSC01 0 

ZS  = 0. 

ZTNSOO 1 1 

C 

ZTNSOO 1 2 

EC  1 1=1, R 

ZTNSOO 1 3 

X F = X ♦ F 

ZTNSOO 14 

XF  P=X*2. »F 

ZTNS0C18 

2E  = ZS*r/J.«  (S  (00,  TO,  P,  X)  *4.*S  (EC.TO.P,  XP)  *S  <E D, TO , R,  XI  E ) ) 

ZTNS0016 

1 X=X*2.*t 

ZTNS001 7 

c 

ZTNSOO 1 8 

EETORN 

ZTNS0019 

END 

7.TNS0020 
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S U RROUT INf  SCOT (ITER, TAC.ANU.T CIS) 
DlflENSICN  XCO?  (5)  ,COF  <c) 
CCNNON/CNF/ROCTH  (4)  , BCCT 1 (4) 
nr*4 

PI  = 3.  14  159 
Xf T=FI*FLCAT (TTER) 

A*  - 4 . / TAC 

F = 6./TA  U«*2-  1./TA'1**4 
C=-4./TAU** I 
CE=XPI/  (2.»ANU*YCIS) 

E= ( 1 .-DD»CD) /T  A U • * 4 

XCOF  (1)  =C 

XCOF  (?)  = C 

XCOF  (3)  -t 

XCOF  ( 4 ) - A 

XCOF  (5)  * 1 . 

CAII.  PCIRT  (XCOF. COP,  FF,  RCOTR. RCCTI. I E R) 

RETURN 

END 


BOOT000 1 
BOOT0002 
BOOTOOO  3 
BOOT0004 
BOOTOOOS 
POOT0006 
HOOT0007 
BOOTOOOS 
ROOT0009 
BOOTOOI 0 
ROOT001 1 
BOOTOOI 2 
BOOTOOI 3 
BOOT  00 1 4 
ROOT0015 
BOOTOO 1 6 
BOOTOOI 7 
BOOTOO 1 8 
BOOTOO 1 9 


C 


c 


PUNi%TIC  N F (D,T,  R,T) 

FNC10001 

A - T ♦ X 

F NC 1 0002 

FNC1310  I 

F f * ? > - a 

FNC10004 

r ! =•*  >S  ( ? F ) 

F NC 1 0005 

E=som  (t  f) 

FNC10006 

A 1 = D*X 

F NC 1 0007 

A 2 = - A 1 

FNC1000H 

IF  ( A2.GT. 50. ) F=0. 

FNC100.79 

IF  (72. GT. 50.)  RETURN 

FNC100 10 

F = R * c X P (11) 

F NC 1 00 1 1 

RETURN 

FNC100 12 

END 

F NC 1 00  1 3 

FI1NCTTCN  S (rt.TO.R.TFTA) 

FNC10014 

A*7ETAMC 

FNC10015 

A R G = A F S ( R • R-  A • A) 

F NC 1 00 1 G 

•R*  IS  IFSS  TUAN  OR  FCAI 

TO  'A*.  PUT  DUE 

TO  THE  SUGHT 

F NC 10017 

NU4F7FC  At  ERROR,  R * R - A • A 

SONFTIfPS  GIVES 

VALUES,  THUS 

F NC 1 00 1 H 

TAKING  TIE  A PS  (4  • R- A * A) 

IS  PREFERRED. 

F NCI oo  in 

11*  SORT  (ARC.) 

F NC 1 0020 

AA.*CD»2ETf 

FNC  10021 

P NC  1 0022 

IF  (AA.LT.-50.)  S = 0. 

FNC10023 

IF  (AA.LT.-5C.)  RETURN 

FNC10024 

FNC 1 0025 

S=  M *7.  F T A * F X P (AA) 

FNC10026 

F FTURN 

FNC 10027 

END 

PNC10020 

I 
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P 0 NCT IC  N PS  (X) 

A = SI  N (X  ) 

E-  X*COS  (X) 

PS=A-S 

fp  turn 
IKD 

PUNCTICH  FC(X) 

» = COS  (X) 

E=X«ST  >1  (X) 

FC=A»E 

RETURN 

END 

P U NCT  IC  N SS  (P  PS. TO,  RC) 
PI=3. 1415928 
IP(RO.fg.C.)  SS=0. 
If(RO.?C-0.)  RPTI1RN 
AlPA-ASIN (EPS/PC) 

E = RO*COS ( / L P f ) 

EE=B/2. 

A 1 = 2.*EFS» (TC-f  F) 

A2- RO**  2* (PI-ALPA) 

SS=A1*»2 

RETURN 

END 

EUNCTTCN  ESIN(X) 

A = SQRT (1.-X»X) 

i = V> 

ASIN=»TAN ( Q) 

RETURN 

END 


FNC2000 1 
FNC20002 
PWC20003 
FNC20004 
PNC20005 
FNC20006 
FNC2000  7 
FNC2000R 
FNC20009 
FNC20U10 
TUC200 1 1 
FNC2001 2 
FNC20C 1 3 
FNC20014 
FNC2001S 
FHC20016 
FNC20017 
FNC20018 
FNC20019 
E NC20020 
FNC2002  t 
FNC20022 
FNC20023 
PHC20024 
FNC20D25 
FNC20028 
PNC  20027 
FNC20028 
FNC20029 
FNC20030 
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